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With numerical experiments we explore the feasibility of using high frequency waves for probing the magnetic fields in 
the photosphere and the chromosphere of the Sun. We track a plane-parallel, monochromatic wave that propagates through 
a non- stationary, realistic atmosphere, from the convection-zone through the photosphere into the magnetically dominated 
chromosphere, where it gets refracted and reflected. We compare the wave travel time between two fixed geometrical 
height levels in the atmosphere (representing the formation height of two spectral lines) with the topography of the surface 
of equal magnetic and thermal energy density (the magnetic canopy or (3 — 1 contour) and find good correspondence 
between the two. We conclude that high frequency waves indeed bear information on the topography of the 'magnetic 
canopy'. 
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1 What is CO s BOLD? 

CO s BOLD is a computer code designed for simulating hy- 
drodynamic processes including radiative transfer in the ou- 
ter and inner layers of stars. Additionally, it can treat the in- 
teraction of magnetic fields with a plasma in the magnetohy- 
drodynamic (MHD) approximation, non-equilibrium chem- 
ical reaction networks, dynamic hydrogen ionization, and 
dust formation in stellar atmospheres. C0 5 BOLD stands for 
Conservative COde for the Computation of COmpressible 
COnvection in a BOx of L Dimensions with L=2,3. 

In the past few years, C0 5 BOLD has been extensively 
used for 'box in a star' computations (for a review see Stef- 
fen 2007), where a small piece of a star's interior and/or 
atmosphere is simulated. A typical result of these kinds of 
computations is the emergent intensity from the stellar sur- 
face. In this way, the surface granulation pattern of the Sun 
and other stars, e.g., A-stars (Kochukhov et al. 2006) or red 
giants have been simulated. In the case of the Sun, a de- 
tailed comparison between high resolution observation and 
simulation is possible and allows for a validation of the lat- 
ter (Steffen, Ludwig, & Freytag 2002). On the other hand, 
CO s BOLD is also used for 'star in a box' computations in 
which the computational box contains an entire star, e.g., the 
red supergiant Betelgeuse (Freytag, Steffen, & Dorch 2002; 
Chiavassa et al. 2006). 

Recent applications of C0 5 BOLD focus on the dynam- 
ical evolution of carbon monoxide in the solar photosphere 



Corresponding author: e-mail: steiner@kis.uni-freiburg.de 



(Wedemeyer-Bohm et al. 2005; Wedemeyer-Bohm & Stef- 
fen 2007) and on the dynamical hydrogen ionization (Leen- 
aarts & Wedemeyer-Bohm 2006), where deviations from 
statistical ionization equilibrium occur. These works have 
been carried out in connection with efforts for a simulation 
including the layers from the convection zone to the chro- 
mosphere of the Sun (Wedemeyer et al. 2004), which now 
have been expanded in order to include the effects of mag- 
netic fields (Schaffenberger et al. 2005, 2006). Details of 
these MHD-simulations are given in the next section. 

CO s BOLD works with Cartesian, non-equidistant but 
static grids, realistic equations of state, multidimensional 
radiation transfer, realistic opacities, and it features vari- 
ous boundary conditions. As a simplest option the radiative 
transfer can be treated frequency-independently ('grey'). 
Alternatively, the frequency-dependence can be taken into 
account by using a multi-group scheme that solves the ra- 
diative transfer problem for different opacity groups, which 
takes effects of line blanketing into account. 

The multidimensional problem is reduced to 1-D prob- 
lems by dimensional splitting. Each of these 1-D problems 
is solved with a Godunov-type finite-volume scheme us- 
ing an approximate Riemann solver modified for a general 
equation of state and gravity. A Roe-type approximate Rie- 
mann solver is used for integration of the hydrodynamic 
equations, whereas in the case of the MHD-equations we 
use a HLL-solver of second order accuracy, instead. More 
details on the MHD-solver can be found in Schaffenberger 
et al. (2005). A good introduction to the above mentioned 
numerical techniques can be found in Laney (1998), LeV- 
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Fig. 1 Left: Magnetic flux density on a color scale, where red signifies strong (greater or approximately 0. 1 T) and blue 
weak fields. The surface of equipartition between magnetic and thermal energy density (]3 = 1) is indicated as blue surface 
(on the lateral boundaries by a yellow curve). The black curve indicates optical depth unity. Right: The temperature in 
colors shows granules and cool intergranular downflows at a height that corresponds to the mean optical depth unity. Note 
the mesh-work like pattern of hot shock waves in the chromospheric layers. From Wedemeyer-Bohm (2007). 



eque (2002), LeVeque et al. (1998), or Toro (1999). CO 5 - 
BOLD is programmed in FORTRAN 90 with OpenMP di- 
rectives. The manual for C0 5 BOLD can be found under 
www.astro.uu.serbf/co5bold_main.html while talks de- 
livered at the first CO s BOLD workshop are at http://folk 
uio.no/svenwe/cws2006/cws_cont.html. For the analysis 
and visualization a CO s BOLD analysis tool (CAT, http:// 
folk.uio.no/svenwe/co5bold/cat.html) exists. 

2 MHD-simulation from the convection zone 
to the chromosphere 

A three-dimensional magnetohydrodynamic simulation of 
the integral layers from the upper convection zone to the 
middle chromosphere was carried out by Schaffenberger et 
al. (2005, 2006). The computational domain in this simu- 
lation extends over a height range of 2800 km, of which 
1400 km reach below the mean surface of optical depth 
unity and 1400 km above it. The horizontal dimensions are 
4800 km x 4800 km. With 120 3 grid cells, the spatial reso- 
lution in the horizontal direction is 40 km, while in the ver- 
tical direction it is 50 km through the convection-zone layer 
and 20 km throughout the photosphere and chromosphere. 
The lateral boundary conditions are periodic in all variables, 
whereas the lower boundary is 'open' in the sense that the 
fluid can freely flow in and out of the computational do- 
main under the condition of vanishing total mass flux. The 
specific entropy of the inflowing mass is fixed to a value pre- 
viously determined so as to yield solar radiative flux at the 
upper boundary. The upper boundary can be chosen to be 
transmitting so that acoustic waves can leave the computa- 
tional domain with no significant reflection at the boundary. 



Stress-free conditions are in effect for the horizontal veloc- 
ities, viz., dv Xt y/dz = 0. 

The MHD simulation starts with a homogeneous, ver- 
tical, unipolar magnetic field of a flux density of 0.001 T 
(10 G) superposed on a previously computed, relaxed model 
of thermal convection. This flux density is thought to mimic 
magnetoconvection in a very quiet network-cell interior. 
The magnetic field is constrained to have vanishing hori- 
zontal components at the top and bottom boundary but lines 
of force can freely move in the horizontal direction, allow- 
ing for flux concentrations to extend right to the boundaries. 
Although this condition is quite stringent, especially at the 
top boundary, it still allows the magnetic field to freely ex- 
pand with height through the photospheric layers into the 
more or less homogeneous chromospheric field, different 
from conventional simulations that extend to a height of typ- 
ically 600 km, only. 

Subsequent to superposition of the magnetic field, flux 
expulsion from the granule centers takes place and within 
less than 5 minutes, the magnetic field concentrates in nar- 
row sheets and small knots near the surface of optical depth 
unity with field strengths up to approximately 1 kG. Occa- 
sionally, these magnetic flux concentrations extend down to 
the bottom boundary at a depth of 1400 km but more often, 
they disperse again at a depth of less than 1000 km leaving 
flux concentrations of a strength of a few hundred Gauss 
only. 

Fig.[T]gives a synoptic summary of phenomena that oc- 
cur in this simulation. The box on the left hand side ren- 
ders the logarithm of the magnetic flux density on a color 
scale, where red signifies strong (greater or approximately 
0.1 T) and blue weak fields. The surface of equipartition 
between magnetic and thermal energy density (f3 = 1) is 
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Fig. 2 Still images from a time series for the instant t = 1368 s after starting with an initial homogeneous vertical field 
of 100 G. Left: Logarithmic magnetic flux density. Right: Logarithm of the ratio of thermal to magnetic energy density 
(plasma-/?). Also shown are the contour of (3 = 1 (thin curve) and the surface of optical depth unity, T500 nm = 1- A strong 
magnetic flux sheet has formed near x = 4000 km that causes the funnel of low (3 visible in the right hand panel. 



indicated as blue surface and on the front side by a yel- 
low curve. It is corrugated and continuously changes shape 
on a time-scale of a few minutes. The black curve indi- 
cates optical depth unity. One can see that close to the front 
side boundary, a magnetic flux concentration has formed. It 
leads to a small dip in the surface of optical depth unity due 
to the reduced density within the flux concentration com- 
pared to the density in the environment at the same geo- 
metrical height. The box on the right hand side shows the 
temperature in colors. We see granules (yellow to light red) 
near the height that corresponds to optical depth unity and a 
network of cool (dark red) intergranular downflows. In the 
chromosphere exists a dynamic, mesh-work like pattern of 
hot shock waves with cool post-shock regions in between. 

3 Local helioseismic experiments with high 
frequency waves 

Classical helioseismic tools usually rely on the analysis of 
waves with frequencies below the acoustic cutoff that are es- 
sentially trapped in the internal acoustic cavity of the Sun. 
These waves are largely evanescent in the outer atmosphere 
and carry little information on this region. The situation 
changes for waves with frequencies above the acoustic cut- 
off frequency that can freely propagate in the atmosphere. 
When they encounter a change in the dispersive character- 
istics of the medium, they are refracted and reflected (Fin- 
sterle et al. 2004). 

This effect has been employed by Finsterle et al. (2004) 
in order to obtain the three-dimensional topography of the 
'magnetic canopy' in and around active regions by deter- 
mining the propagation behavior of high-frequency acoustic 
waves in the solar chromosphere. Their method bears con- 
siderable potential for the exploration of the magnetic field 
in the solar atmosphere by helioseismological means. 



In order to test the reliability and to further explore and 
assess this method we have carried out a first series of two- 
dimensional MHD-simulations of the propagation of high 
frequency waves through a magnetically structured, realis- 
tic atmosphere. We start with a simulation very similar to 
that described in the previous section, but carried out in a 
two-dimensional computational domain of 4900 km width 
and spanning a height range of 2900 km, of which 1300 km 
reach into the convection zone and 1600 km above the aver- 
age level of optical depth unity. 

Since we start with a mean magnetic flux density of 
10 G and 100 G, the magnetic structures that form in the 
course of time represent magnetic fields in a very quiet net- 
work cell interior and in the magnetic network, respectively, 
rather than the active region fields observed by Finsterle et 
al. (2004). Also, different from active region fields, this field 
evolves on the much shorter time scale of granular evolu- 
tion. Therefore, and for reasons of computational costs, we 
do not analyze a long duration time series as is normally 
necessary in local helioseismology. Instead, we introduce 
a velocity perturbation of given frequency and amplitude at 
the bottom of the computational domain with which we gen- 
erate a monochromatic, plane parallel wave that propagates 
within a time span of about 200 s across the height range of 
the computational domain. Experiments with various ampli- 
tudes and frequencies have been carried out. In the follow- 
ing, however, we only use the values stated in connection 
with Eq. ([]])■ 

We then determine the arrival time of the wave front at 
three different heights in the atmosphere that should roughly 
represent the heights of maximal Doppler response of the 
spectral lines Ni I 676.8 nm at 200 km, K I 769.9 nm at 
420 km, and Na 1 589.0 nm at 800 km (Finsterle et al. 2004). 
Subsequently, we compute the travel times between these 
heights. They are thought to become modified through the 
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presence of a strong magnetic field in the sense that a strong 
magnetic field reduces the travel time. At this stage we do 
not carry out radiation transfer for computing the response 
of spectral lines to the perturbation but simply determine the 
velocity at the corresponding (fixed) line-formation heights 
in the computational box for determining the arrival time of 
the wave front. The wave front is detected when the velocity 
perturbation surpasses a height dependent velocity thresh- 
old of less than 100 m/s. 

In the following we consider the time span from t = 
1200 s to t = 1400 s of a simulation that started with a 
homogeneous vertical field of 100 G at t = s. Fig. 
shows a snapshot of the instant t = 1368 s. The panel on 
the left hand side shows the logarithmic magnetic flux den- 
sity, where the values in Gauss are indicated in the grey- 
scale bar in the top margin of the figure. A strong magnetic 
flux sheet has formed near x = 4000 km. It leads to a dip 
in the surface of optical depth unity, visible in the panel on 
the right hand side. There, also the contour of {3 = 1, where 
(3 = Pgas/ {B 2 /8tt), i.e., the contour of equipartition of ther- 
mal and magnetic energy density, is indicated. The magnetic 
flux concentration causes a funnel of low (3. 

The (3=1 contour delimits the transition from the weak 
to the strong field regime. In the latter regime the magnetic 
field dominates over thermal pressure. The gas pressure ex- 
ponentially decreases with height in a gravitationally strat- 
ified atmosphere but the magnetic pressure can only de- 
crease to a lower limit set by the initial average flux den- 
sity. Once the magnetic field, which may be concentrated 
in small tubes and sheets in the deep photosphere, has ex- 
panded to completely fill the available space, it cannot fur- 
ther expand and its strength remains essentially constant 
with height. Therefore, (3 unabatedly decreases with height 
so that with increasing height the atmosphere becomes 
sooner or later magnetically dominated. Correspondingly, 
the (3 = 1 contour in Fig. [2] (right) extends in a more or less 
horizontal direction at a height of 500 km with the excep- 
tions of the location of the strong flux sheet, where it dips 
even below z = 0, and two islands higher up in the atmo- 
sphere, caused by chromospheric shock waves. 

The (3 = 1 contour also marks the region of wave trans- 
mission and conversion (Rosenthal et al. 2002; Bogdan et 
al. 2003; Cally 2007). For example an incident acoustic 
plane wave converts to a fast magnetic wave and gets re- 
flected back into the atmosphere. This effect is visible in 
Fig- El which shows, from bottom to top, a time sequence of 
a plane wave traveling through the magnetically structured 
and simultaneously evolving atmosphere of Fig. [2] In order 
to visualize the wave we have run the simulation twice from 
t = 1200 s to t = 1400 s, once without perturbation and 
once with a velocity perturbation at the bottom of the com- 
putational domain of the form 

v z (t) = w sin(27r(i - t )v) , (1) 

where vo was chosen to be 0.05 km/s and v = 20 mHz. 
The two runs were carried out with identical time stepping. 
Following that, the two velocity fields are subtracted, which 
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Fig. 3 A plane parallel wave with frequency 20 mHz trav- 
els through convecting plasma into the magnetically struc- 
tured photosphere and further into the low (3 (magnetically 
dominated) chromosphere. The three panels show the dif- 
ference in absolute velocity between the perturbed and the 
unperturbed solution 100 s, 148 s, and 168 s (from bottom 
to top) after the start of the perturbation (Eq. Q]). Magnetic 
field and plasma (3 corresponding to the instant of the top 
panel are given in Fig. [2] A spurious velocity signal that is 
believed to have numerical origin dominates the very low (3 
region in the upper right and upper left corners. The hori- 
zontally running black curve near z = indicates optical 
depth T5oonm = 1- The velocity scaling is logarithmic with 
dimension [cm/s]. 
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then reveals the wave perturbation that travels on top of the 
non-stationary evolving convective motion. 

Due to a numerical problem of unknown origin, consid- 
erable spurious velocities occur in the region of very low 
plasma- (3 in the chromospheric layers of the magnetic flux 
concentration. This problem will be addressed in a next sta- 
ge of the present project. Nevertheless, we can clearly ob- 
serve the following phenomena. Initially the wave front is 
retarded at the location of the flux concentration (at x « 
4000 km) because of the lower temperature of the plasma 
within the flux concentration, hence the lower sound speed 
(bottom panel, 100 s after start of the perturbation). It un- 
dergoes acceleration when entering the funnel of (3 < 1 so 
that the wave within the funnel is now leading 148 s after 
start of the perturbation (middle panel). At the same time 
the wave becomes fast magnetic in character and starts to 
refract. Thus, the wave front becomes inclined until aligned 
with the vertical direction and it further turns until trav- 
eling into the downward direction. In the top panel, only 
20 s later, the wave front in the low-/? funnel (extending 
from (x, z) = (2400, 1500) to (3400, 500)) has already 
completely turned around and travels back into the atmo- 
sphere. Therefore, we can speak of a reflection of the mag- 
netic wave in the low-/? region. A similar fanning out of the 
wave front occurs around x = 1400 km when it reaches the 
'canopy' -height where (3 = 1 at about 500 km. 

Clearly, the travel time for the wave in the low-/3 re- 
gion is much smaller than elsewhere. Fig.|4]shows the travel 
time between the line formation heights of Ni I, 676.8 nm at 
200 km and K I 769.9 nm at 420 km (thick solid curve) as 
a function of horizontal distance, x. The time axis is given 
on the left hand side. Superposed on this plot is the contour 
of (3 = 1 (dash-dotted curve), where the height in the atmo- 
sphere is given by the z-axis on the right hand side. Clearly, 
the travel-time curve follows the dip of the /^-contour at the 
location of the flux concentration, which demonstrates that 
a mapping of the (3=1 surface by helioseismic methods is 
in principle possible. 

Because of the spurious velocity noise in the region of 
very low (3, we were unable to reliably determine the wave- 
front arrival at the line-formation height of Na I D2, 589.0 
nm at 800 km. But we expect for the travel time difference 
between 800 and 420 km an even more pronounced dip at 
the location of the magnetic flux concentration than for the 
difference between 200 and 420 km. 

We have also carried out a similar simulation with an 
initial vertical field of 10 G flux density. Because of the 
lower field strength the average height where (3=1 (canopy 
height) is in this case at approximately 1000 km compared 
to 500 km in the case of 100 G. We could think of the two 
simulations with 10 G and 100 G initial field strength rep- 
resenting very quiet network interior and network regions, 
respectively. Examination of the travel times for 20 mHz 
waves reveals no major difference in the layer from 200 km 
to 420 km. In the 100 G case the mean travel time is 23.9 s 
and 25.7 s when excluding the strong flux concentration that 
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Fig. 4 Wave travel time across the layer from z = 200 km 
to z = 420 km as a function of horizontal distance (thick 
solid curve). Superposed is the contour of (3 = 1 (magnetic 
and thermal equipartition), for which the height is indicated 
in the right hand side ordinate (dash-dotted curve). Note that 
the travel time markedly decreases where the low (3 region 
intrudes this layer. 



has formed at x rj 4000 km. The corresponding time in the 
10 G case is 28.0 s. Since this layer is below the canopy in 
both cases, there is only a small difference in wave travel- 
time. However, for the layer from 420 km to 800 km we 
obtain in the 100 G case a travel time of 25.1 s (when ex- 
cluding the strong flux concentration, where reliable values 
are unavailable). The corresponding travel time in the 10 G 
case is 39.2 s, much longer because in this case this layer is 
still below the low-/3 regime while it is entirely inside it in 
the 100 G case. Again this demonstrates that measurements 
of the wave travel-time can differentiate between network 
and internetwork regions. 

4 Conclusions 

We have carried out local helioseismic experiments using 
the magnetohydrodynamic code C0 5 BOLD. A plane-paral- 
lel monochromatic wave perturbation that is generated at 
the bottom of the computational domain within the convec- 
tion zone propagates into a non-stationary, realistic atmo- 
sphere that extends up into the middle chromosphere. There, 
it enters a regime in which the magnetic energy density sur- 
passes the thermal energy density and converts from acous- 
tic to magnetic in nature. It is experimentally demonstrated 
how a magnetic flux concentration leads to refraction and 
reflection of the wave. 

The experiments also reveal a numerical problem of un- 
known origin. Two runs that are identical with the exception 
of the wave perturbation at the bottom of the computational 
domain show considerable differences in the velocity field 
in regions of (3 « 10~ 2 from the very beginning. In these 
tenuous regions, smallest pressure imbalances that may oc- 
cur because of the thermal energy density being only a small 
fraction of the total energy density there, must lead to large 
accelerations. 

We show that the wave travel-time between two fixed 
levels in the atmosphere bears information on the nature 
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of the wave and consequently on the magnetic field. The 
travel time is reduced in regions of low (3 (strong magnetic 
field). For a monochromatic wave of frequency 20 mHz we 
demonstrate with Fig. @] that the travel time between the 
heights of 200 and 420 km in seconds times 15 matches 
approximately the height of the (3=1 surface in km. In 
particular, the travel-time curve delineates a funnel of low (3 
that is caused by a local magnetic flux concentration. 

The region and magnetic structure considered in the pre- 
sent simulation is much smaller than the active region ob- 
served by Finsterle et al. (2004) and we use 20 mHz waves, 
instead of the 7 mHz employed by Finsterle et al. (2004). 
Also, by using plane parallel waves we assume the wave 
coherence to be always larger than the magnetic structure 
under investigation. Subject to these reservations, the nu- 
merical experiments carried out in this paper support the 
conclusions of Finsterle et al. (2004) and their proposition 
using high frequency waves for mapping the magnetic to- 
pography in the chromosphere. 
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